{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-d485c91a8c052319",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "# Exercise 8: Supervised Learning"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-ef5b12a5b0a9750d",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "In this exercise a few basic supervised learning pipelines are elaborated.\n",
    "* Data analysis\n",
    "* Feature engineering\n",
    "* Cross-validation\n",
    "* Regression\n",
    "  * Linear models and neural networks\n",
    "* Classification\n",
    "  * Linear models and neural networks\n",
    "\n",
    "Moreover, working with the popular __pandas__ library is presented here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-47de89ead89817f8",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.colors import rgb2hex\n",
    "import seaborn as sns\n",
    "from sklearn.preprocessing import StandardScaler, MinMaxScaler\n",
    "plt.style.use('default')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "jp-MarkdownHeadingCollapsed": true
   },
   "source": [
    "## a) Data Set"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-e104b64c287bccc1",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "In this exercise we will use a [data set hosted on Kaggle](https://www.kaggle.com/hankelea/system-identification-of-an-electric-motor).\n",
    "For ease of use and simplicity, we are going to work on a reduced version though, that is provided with this notebook already.\n",
    "\n",
    "The data set represents lab measurements from a test bench, where a three-phase electric motor was externally speed-controlled while being induced with different voltages, leading to corresponding currents and, thus, torque on the motor. In terms of reinforcement learning the electric motor can be understood as the environment, the observed electrical state of the motor as the state vector and the applied voltages/discrete switching states of the converter as the actions.\n",
    "\n",
    "### Context\n",
    "The most important aspect of electric vehicles is their efficiency or achievable range. \n",
    "In order to achieve high efficiency, it is essential to avoid the over-dimensioning of the drive train during the design phase.\n",
    "This can only be achieved if the dynamic behavior of the drive train is accurately known by the controller/agent. \n",
    "The task of the controller is to achieve a desired torque at the wheels of the car by controlling the currents of the electric motor.\n",
    "\n",
    "Control, in this context, is about choosing appropriate, discrete switching states for the transistors in the attached power electronic unit, that affect the voltage applied to the motor.\n",
    "Knowing which switching state leads to what current change at each point in time helps the controller to choose an appropriate action.\n",
    "Since in the real world many parasitic, intangible and volatile effects are imposed on currents formation, modeling them statistically is an auspicious endeavor.\n",
    "\n",
    "### Features\n",
    "\n",
    "Most of the time, a data set is a table/matrix/2d-tensor with labeled columns (and rows).\n",
    "\n",
    "Here, rows are consecutively enumerated while columns are the following:\n",
    "\n",
    "* __n_1k__: elementary vector applied between (k-1) and k. The elementary vector denotes the label-encoded switching state in the converter.\n",
    "* __id_k__: d-current at k in A (Ampere).\n",
    "* __iq_k__: q-current at k in A (Ampere). \n",
    "* __epsilon_k__: rotor angle at k in rad (Radian).\n",
    "* __n_k__: elementary vector applied between k and (k+1). \n",
    "* __id_k1__: d-current at (k+1) in A (Ampere).\n",
    "* __iq_k1__: q-current at (k+1) in A (Ampere).\n",
    "\n",
    "[Read up on d/q-transformation](https://en.wikipedia.org/wiki/Direct-quadrature-zero_transformation) to understand the coordinate notation of the currents (assuming the amplitude-invariant transformation than the power-invariant).\n",
    "\n",
    "The enlisted features will be utilized for either regression or classification.\n",
    "Depending on the task, the features denoting input and target quantities will vary, such that the specific mapping will be discussed later.\n",
    "\n",
    "Lab measurements from a motor are commonly time series data.\n",
    "Yet, here the records are shuffled, and we only have the link between two elementary vectors (switching states/ actions) in a row itself.\n",
    "\n",
    "### Optional Literature\n",
    "If you find the following exercise intriguing and you would like to keep on investigating, have a look at the scientific papers that describe the background to the underlying data set:\n",
    "* [Part 1](https://arxiv.org/pdf/2003.07273.pdf)\n",
    "* [Part 2](https://arxiv.org/pdf/2003.06268.pdf)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-b27844f3da8cecb3",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "### Load CSV data set"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-422cb7f5e6909a9c",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "The pandas library makes reading csv files easy, next to its superb data wrangling utilities. It is also able to read and save archived tabular data as in this exercises' case.\n",
    "\n",
    "Think of pandas as a layer on top of numpy, that is optimized for tabular data and attaches column/row names to a matrix.\n",
    "\n",
    "The main pandas class is the so-called _DataFrame_ which holds a table, and the _Series_  that denotes a single column.\n",
    "Both classes do not only consist of arrays representing the data itself but also of an index (row labeling), that affects basic arithmetic operations between different instantiated objects.\n",
    "\n",
    "Please refer to the pandas notebooks from exercise 1 for a refresher."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-37acd3211bc3591f",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "df = pd.read_csv('reduced_emotor.zip')\n",
    "print(df.shape)\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-7808b756256de207",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "### Adding new features\n",
    "There are two main ways to add further features to the dataframe."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-fb6e9f4fc33dd9a5",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# new feature \"pairs\"\n",
    "df = df.assign(pairs=lambda r: r.n_1k.astype(str)+'->'+r.n_k.astype(str))\n",
    "# alternatively\n",
    "df['pairs_alternative'] = df.n_1k.astype(str) + '->' + df.n_k.astype(str) \n",
    "\n",
    "assert all(df.pairs == df.pairs_alternative)\n",
    "# remove a feature\n",
    "del df['pairs_alternative']\n",
    "# alternatively\n",
    "# df = df.drop(['pairs_alternative'], axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-aec6adffa5403150",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "df.columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-a3eeef12cf6c43c8",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "df.pairs.value_counts().sort_index().head(10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-929235a4f8d308fb",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "We recognize, that all transitions within a row appear uniformly throughout the dataset, which is not the case for the original kaggle data but has been intentionally rectified for this exercise's reduced version."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "jp-MarkdownHeadingCollapsed": true
   },
   "source": [
    "## b) Exploratory Data Analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-2a311570ae5f45b5",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "Before any feature engineering should be conducted, it is advisable to investigate the general structure of the data.\n",
    "\n",
    "This process is called exploratory data analysis (EDA) and usually comprises \n",
    "* distribution visualizations, \n",
    "* missing data detection, \n",
    "* substitution of such, \n",
    "* linear correlation analysis,\n",
    "* time series analysis\n",
    "* etc."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-b3d8f77b0430d827",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "### Distribution visualizations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-f5ad227060c3bfb8",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "def plot_dist(_df):\n",
    "    fig, axes = plt.subplots(1,2,sharex=True, sharey=True)\n",
    "    for c, ax in zip(['n_1k', 'n_k'], axes.flatten()):\n",
    "        sns.countplot(x=c, data=_df,  ax=ax)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-c39f8271e9fcbf72",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "plot_dist(df)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-dbd820d77bc013a6",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "Each elementary vector occurs around 350 thousand times for both, at time _k_ and time _k-1_ .\n",
    "\n",
    "In the following we plot the distribution of the currents and the motor angle unfolded over each elementary vector.\n",
    "\n",
    "For this we utilize pandas' _groupby()_ function, which helps to group (and even aggregate) the complete data set with respect to a categorical feature."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-773db26050061bc0",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "analyzed_cols = [c for c in df if c not in ['n_k', 'n_1k', 'pairs']]\n",
    "\n",
    "# prepare colors\n",
    "color_list = plt.cm.Set3(np.linspace(0, 1, 12))\n",
    "color_d = {'id_k': rgb2hex(color_list[0]),\n",
    "           'id_k1': rgb2hex(color_list[0]),\n",
    "           'iq_k': rgb2hex(color_list[1]),\n",
    "           'iq_k1': rgb2hex(color_list[1]),\n",
    "           'epsilon_k': rgb2hex(color_list[2])}\n",
    "\n",
    "def dist_plot(dframe, cols_to_plot):\n",
    "    unique_elem_vecs = dframe['n_k'].nunique()\n",
    "    # create plot\n",
    "    fig, axes = plt.subplots(nrows=unique_elem_vecs, ncols=len(cols_to_plot), \n",
    "                             sharex='col', sharey='col', figsize=(3*unique_elem_vecs, 10))\n",
    "    for k, _df in dframe.groupby('n_k'):\n",
    "        for i, c in enumerate(cols_to_plot):\n",
    "            sns.histplot(_df[c], ax=axes[k-1, i], color=color_d.get(c, None))\n",
    "            if i == 0:\n",
    "                axes[k-1, i].set_ylabel(f'n_k = {k}')\n",
    "            if k == 1:\n",
    "                axes[k-1, i].set_title(c)\n",
    "            xlbl = c if k == 7 else ''\n",
    "            axes[k-1, i].set_xlabel(xlbl)\n",
    "\n",
    "    plt.tight_layout()\n",
    "dist_plot(df, analyzed_cols)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-0246df5fc7d9b075",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "_(Double-click picture if it appears too small)_.\n",
    "\n",
    "It seems that _epsilon_k_ exhibits a sine shaped histogram while its range covers $[-\\pi, \\pi]$.\n",
    "\n",
    "More subtle, we recognize a semi-sphere shape of the 2d histogram between the currents (remember, d and q currents are to be plotted perpendicular to each other).\n",
    "It might be auspicious, to add another feature denoting the current vector norm sqrt(id^2 + iq^2).\n",
    "From the principle of d/q-transformation, we know that the $\\ell_2$-norm of the d- and q-current denotes the current amplitude in the system.\n",
    "From that we recognize, that smaller current magnitudes seem more likely than higher."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-7368c69fea2082d4",
     "locked": true,
     "points": 0,
     "schema_version": 3,
     "solution": false,
     "task": true
    }
   },
   "source": [
    "#### Task 1: Add features and plot their distributions\n",
    "* Add sine and cosine of the rotor angle to the dataframe\n",
    "* Add the current vector norm of both time steps to the dataframe: $i_{norm} = \\sqrt{i_d^2 + i_q^2}$\n",
    "* Plot their distributions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": true,
     "grade_id": "cell-273354ee2c57bdb1",
     "locked": false,
     "points": 0,
     "schema_version": 3,
     "solution": true,
     "task": false
    }
   },
   "source": [
    "#### Solution 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbgrader": {
     "grade": true,
     "grade_id": "cell-dbd6ba2be898d482",
     "locked": false,
     "points": 0,
     "schema_version": 3,
     "solution": true,
     "task": false
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# YOUR CODE HERE\n",
    "raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbgrader": {
     "grade": true,
     "grade_id": "cell-f7e3783625157297",
     "locked": false,
     "points": 0,
     "schema_version": 3,
     "solution": true,
     "task": false
    },
    "scrolled": true,
    "tags": []
   },
   "outputs": [],
   "source": [
    "dist_plot(df, [c for c in df if c not in ['n_k', 'n_1k', 'pairs']])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-2e6ec9eb45edb784",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "### Scatter plots\n",
    "Probably, the previous currents are good indicators for the next currents.\n",
    "Let's visualize this relationship unfold over each transition."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-d3b911cceb0a7f66",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    },
    "scrolled": true,
    "tags": []
   },
   "outputs": [],
   "source": [
    "n_cols = 7\n",
    "n_rows = 7\n",
    "fig, axes = plt.subplots(nrows=n_rows, ncols=n_cols, \n",
    "                         sharex=True, sharey=True, figsize=(20, 20))\n",
    "df = df.sort_values(by='pairs')\n",
    "for ax, (k, _df) in zip(axes.flatten(), df.groupby('pairs')):\n",
    "    t_prev, t = k.split('->')\n",
    "    ax.scatter(_df.id_k, _df.id_k1, s=.5, label='d')\n",
    "    ax.scatter(_df.iq_k, _df.iq_k1, s=.5, label='q')\n",
    "    if t == '1':\n",
    "        ax.set_ylabel(f'n_1k = {t_prev}\\n i_k1 in A')\n",
    "    if t_prev == '7':\n",
    "        ax.set_xlabel(f'i_k in A\\n n_k = {t}')\n",
    "    ax.grid(alpha=.5)\n",
    "plt.tight_layout()\n",
    "\n",
    "ax.legend(ncol=2, loc='lower center', bbox_to_anchor=(.5, 1), bbox_transform=fig.transFigure, markerscale=6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-d2c046101730d429",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "Apparently, transitions to switching_state == 1 come with far less scatter than for other switching states.\n",
    "Beyond this, the overall scatter seems relatively linear.\n",
    "\n",
    "We investigate this by a linear correlation heat map."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-845c3e1862504d6a",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "### Linear Correlation\n",
    "The [Pearson correlation coefficient](https://en.wikipedia.org/wiki/Pearson_correlation_coefficient) indicates the strength of the linear correlation between two random variables."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-e6507791679213e1",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "corr = df.drop(['pairs'], axis=1).corr()\n",
    "# Generate a mask for the upper triangle\n",
    "mask = np.zeros_like(corr, dtype=bool)\n",
    "mask[np.triu_indices_from(mask)] = True\n",
    "\n",
    "plt.figure(figsize=(8,8))\n",
    "_ = sns.heatmap(corr, mask=mask, center=0, square=True, linewidths=.5, cbar_kws={\"shrink\": .8})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-45437311f7e70024",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "It becomes obvious, that features from previous time steps exhibit the strongest linear correlation with those of the actual time step.\n",
    "\n",
    "This suggests a strong timely dependency, currents follow a [non-stationary process](https://en.wikipedia.org/wiki/Stationary_process).\n",
    "\n",
    "It is very likely, that more history information for each point in time would reveal precious patterns for an exact forecast of the currents, but these are not available here."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "jp-MarkdownHeadingCollapsed": true,
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-3d387b56880b066c",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "## c) Cross Validation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-dde17289924955d7",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "Before we start modeling, we need to agree on a certain cross-validation (CV) technique.\n",
    "We will choose a method from scikit-learn.\n",
    "\n",
    "Please install it with\n",
    "```\n",
    "pip install scikit-learn\n",
    "```\n",
    "\n",
    "Since distribution of elementary vectors is rather balanced, we will go with a simple CV: 5-Fold CV.\n",
    "\n",
    "Recall from the lecture what K-Fold CV stands for."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-d3c24d634e4afe84",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "![](./kfold-cv.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-337cbc5aaf52aad6",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "from sklearn.model_selection import KFold"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-c067bf5e788d3108",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "kf = KFold(n_splits=5, shuffle=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "jp-MarkdownHeadingCollapsed": true
   },
   "source": [
    "## d) Regression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-613a9001732e0616",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "We start off with a regression example: Given the data set, predict the future currents at any time.\n",
    "\n",
    "### Linear Regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-3022be79207dcf45",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "from sklearn.linear_model import LinearRegression\n",
    "from sklearn.metrics import mean_squared_error as mse\n",
    "\n",
    "\n",
    "# more feature engineering\n",
    "enriched_df = df.assign(# The following are non-task-specific featurs\n",
    "                        id_k_sqrd=lambda x: x.id_k**2,\n",
    "                        iq_k_sqrd=lambda x: x.iq_k**2,\n",
    "                        iq_x_id_k=lambda x: x.id_k*x.iq_k,\n",
    "                        i_norm_sqrd=lambda x: x.i_norm_k**2,\n",
    "                        id_k_log=lambda x: np.log1p(x.id_k.abs()),\n",
    "                        iq_k_log=lambda x: np.log1p(x.iq_k),\n",
    "                        id_k_exp=lambda x: np.exp(-x.id_k.abs()),\n",
    "                        iq_k_exp=lambda x: np.exp(-x.iq_k.abs()),\n",
    "                        i_norm_log=lambda x: np.log1p(x.i_norm_k),\n",
    "                        i_norm_exp=lambda x: np.exp(x.i_norm_k.abs()),\n",
    "                        id_plus_iq=lambda x: x.id_k+x.iq_k,\n",
    "                        id_minus_iq=lambda x: x.id_k-x.iq_k,\n",
    "                        id_over_iq=lambda x: x.id_k/ (x.iq_k + 1),\n",
    "                        iq_over_id=lambda x: x.iq_k/ (x.id_k - 1),\n",
    "    \n",
    "                        # The following features encompass the whole data set. Is this problematic?\n",
    "                        #id_pairgroup_normed=lambda x: x.id_k - df.groupby('pairs')['id_k'].transform('mean'),\n",
    "                        #iq_pairgroup_normed=lambda x: x.iq_k - df.groupby('pairs')['iq_k'].transform('mean'),\n",
    "                        #id_ngroup_normed=lambda x: x.id_k - df.groupby('n_k')['id_k'].transform('mean'),\n",
    "                        #iq_ngroup_normed=lambda x: x.iq_k - df.groupby('n_k')['iq_k'].transform('mean'),\n",
    "    \n",
    "                        # The following are features taken from https://ieeexplore.ieee.org/document/9545442\n",
    "                        # You can try them for your neural network training. They should significantly boost performance\n",
    "                        id_k_cos=lambda x : x.id_k * x.cos_eps_k,\n",
    "                        iq_k_cos=lambda x : x.iq_k * x.cos_eps_k,\n",
    "                        id_k_sin=lambda x : x.id_k * x.sin_eps_k,\n",
    "                        iq_k_sin=lambda x : x.iq_k * x.sin_eps_k,\n",
    "                        id_k_sqrd_sin=lambda x : x.id_k_sqrd * x.sin_eps_k,\n",
    "                        iq_k_sqrd_sin=lambda x : x.iq_k_sqrd * x.sin_eps_k,\n",
    "                        id_k_sqrd_cos=lambda x : x.id_k_sqrd * x.cos_eps_k,\n",
    "                        iq_k_sqrd_cos=lambda x : x.iq_k_sqrd * x.cos_eps_k,\n",
    "                        iq_x_id_k_cos=lambda x : x.iq_x_id_k * x.cos_eps_k,\n",
    "                        iq_x_id_k_sin=lambda x : x.iq_x_id_k * x.sin_eps_k,\n",
    "                        sin_2eps_k=lambda x : np.sin(2*df.epsilon_k), \n",
    "                        cos_2eps_k=lambda x : np.cos(2*df.epsilon_k), \n",
    "                        id_k_sin_2eps_k=lambda x : x.id_k * x.sin_2eps_k, \n",
    "                        id_k_cos_2eps_k=lambda x : x.id_k * x.cos_2eps_k, \n",
    "                        iq_k_sin_2eps_k=lambda x : x.iq_k * x.sin_2eps_k, \n",
    "                        iq_k_cos_2eps_k=lambda x : x.iq_k * x.cos_2eps_k, \n",
    "                      )\n",
    "\n",
    "y_cols = ['id_k1', 'iq_k1', 'i_norm_k1']\n",
    "x_cols = [c for c in enriched_df if c not in y_cols + ['pairs'] ]\n",
    "\n",
    "for i, (train_index, test_index) in enumerate(kf.split(enriched_df)):\n",
    "    x_train, y_train = enriched_df.loc[train_index, x_cols], enriched_df.loc[train_index, y_cols]\n",
    "    x_test, y_test = enriched_df.loc[test_index, x_cols], enriched_df.loc[test_index, y_cols].values\n",
    "    \n",
    "    # feature aggregations over the training set\n",
    "    model = LinearRegression()\n",
    "\n",
    "    model.fit(x_train, y_train)\n",
    "    prediction = model.predict(x_test)\n",
    "\n",
    "    result_string = ', '.join(f'{y_cols[j]}: {mse(y_test[:, j], prediction[:, j]):.2f} A²' for j in range(len(y_cols)))\n",
    "    print(f'Fold {i}: \\nMSE:\\t', result_string)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-6e815061baed0ec3",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "The error is homogenous across folds, which affirms our CV strategy.\n",
    "\n",
    "However, the estimation error is fairly large."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-e2446a51101d2445",
     "locked": true,
     "points": 0,
     "schema_version": 3,
     "solution": false,
     "task": true
    }
   },
   "source": [
    "#### Task 2: What did we miss to do before fitting the linear model?\n",
    "Two important preprocessing steps were missed before conducting model training.\n",
    "Which might they be? It has something to do with categorical values and value ranges.\n",
    "\n",
    "In the following, rectify this blunder, and prepare the data appropriately.\n",
    "* Conduct the CV again.\n",
    "* Plot the estimation error of the last fold over the corresponding ground truth.\n",
    "\n",
    "Can you think of more features? Do you think there might be a problem with those features that were aggregated over the full data set in the beginning?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": true,
     "grade_id": "cell-bca37bee3beb1da5",
     "locked": false,
     "points": 0,
     "schema_version": 3,
     "solution": true,
     "task": false
    }
   },
   "source": [
    "#### Solution 2:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbgrader": {
     "grade": true,
     "grade_id": "cell-d9f9c645f73de47c",
     "locked": false,
     "points": 0,
     "schema_version": 3,
     "solution": true,
     "task": false
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Feel free to use the following code for your solution. It is needed for subsequent cells to work\n",
    "# one-hot-encode categoricals\n",
    "ohe_df = enriched_df.assign(**{**{f'n_k_{i}': (enriched_df.n_k == i).astype(int) for i in range(1, 8)},\n",
    "                               **{f'n_1k_{i}': (enriched_df.n_1k == i).astype(int) for i in range(1, 8)}})\\\n",
    "                        .drop(['n_k', 'n_1k', 'pairs'], axis=1)\n",
    "\n",
    "\n",
    "\n",
    "y_cols = ['id_k1', 'iq_k1', 'i_norm_k1']\n",
    "x_cols = [c for c in ohe_df if c not in y_cols]\n",
    "\n",
    "\n",
    "x_cols_to_scale = [c for c in x_cols if not c.startswith(('n_k', 'n_1k'))]\n",
    "\n",
    "# YOUR CODE HERE\n",
    "raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": true,
     "grade_id": "cell-c141066ee4b61c0a",
     "locked": false,
     "points": 0,
     "schema_version": 3,
     "solution": true,
     "task": false
    }
   },
   "source": [
    "### Residuals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbgrader": {
     "grade": true,
     "grade_id": "cell-100b0c82137fbaa7",
     "locked": false,
     "points": 0,
     "schema_version": 3,
     "solution": true,
     "task": false
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# residuals\n",
    "fig, axes = plt.subplots(nrows=1, ncols=3, figsize=(10, 3), sharey='row')\n",
    "error = y_test - prediction\n",
    "\n",
    "for i in range(len(y_cols)):\n",
    "    axes[i].scatter(y_test[:, i], error[:, i], s=.5, label=y_cols[i])\n",
    "    axes[i].grid(alpha=.5)\n",
    "    axes[i].axhline(ls='--')\n",
    "    if i == 0:\n",
    "        axes[i].set_ylabel('Error in A')\n",
    "    axes[i].set_xlabel(y_cols[i]+' in A')\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": true,
     "grade_id": "cell-aa8d7b637f1fee4b",
     "locked": false,
     "points": 0,
     "schema_version": 3,
     "solution": true,
     "task": false
    }
   },
   "source": [
    "The residuals show strong geometrical shapes along the ground truth. Hence, systematic errors are still included."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## e) Neural Networks"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-c095dbc26d1cee80",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "In the following, we showcase how a neural network can be trained to solve the aforementioned prediction problem.\n",
    "\n",
    "Is a neural network better?\n",
    "\n",
    "Be aware that training time of a neural network is substantially longer than for linear models. If you face out-of-memory errors, please reduce the n_units parameter, the batch_size parameter, or subsample the data set.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-2279090db77d2faa",
     "locked": true,
     "points": 0,
     "schema_version": 3,
     "solution": false,
     "task": true
    }
   },
   "source": [
    "#### Bonus task: Try to tweak the hyperparameters of the NN in order to achieve a better performance.\n",
    "\n",
    "If you find a decent setting, you are invited to post your solution to [kaggle kernels](https://www.kaggle.com/hankelea/system-identification-of-an-electric-motor/kernels). Don't be shy!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import torch\n",
    "import torch.nn as nn\n",
    "import torch.optim as optim\n",
    "from torch.utils.data import DataLoader, TensorDataset\n",
    "from tqdm.notebook import tqdm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class FeedForwardNetwork(nn.Module):\n",
    "    \"\"\"Example definition for a very basic feed forward neural network.\"\"\"\n",
    "\n",
    "    def __init__(self, input_dim):\n",
    "        \"\"\"Initialization function for the network.\n",
    "        \n",
    "        This function stacks multiple modules back to back to form a layer stack. The \n",
    "        output of each module is used as the input to the next module.\n",
    "        \"\"\"\n",
    "        super(FeedForwardNetwork, self).__init__()\n",
    "        self.layer_stack = nn.Sequential(\n",
    "            nn.Linear(input_dim, 32),\n",
    "            nn.ReLU(),\n",
    "            nn.Linear(32, 3)\n",
    "        )\n",
    "\n",
    "    def forward(self, x):\n",
    "        \"\"\"Forward pass through the network.\n",
    "\n",
    "        Two things happen here simultaneously:\n",
    "         1. The input x is passed through all modules in the layer stack to compute the\n",
    "             output of the network.\n",
    "         2. The operations gradient functions are stored in a computational graph to\n",
    "             compute the gradient of output w.r.t. the parameters of the network in\n",
    "             the backward pass.\n",
    "    \n",
    "        https://pytorch.org/tutorials/beginner/basics/autogradqs_tutorial.html#more-on-computational-graphs\n",
    "        \"\"\"\n",
    "        return self.layer_stack(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def prepare_data(train_set, test_set, y_cols, x_cols, x_cols_to_scale):\n",
    "    \"\"\"Scales the data and converts it to torch tensors.\"\"\"\n",
    "    \n",
    "    x_scaler, y_scaler = MinMaxScaler(), MinMaxScaler()\n",
    "    x_train = train_set.loc[:, x_cols]\n",
    "    x_train.loc[:, x_cols_to_scale] = x_scaler.fit_transform(x_train.loc[:, x_cols_to_scale])\n",
    "\n",
    "    y_train = y_scaler.fit_transform(train_set.loc[:, y_cols])\n",
    "    x_test = test_set.loc[:, x_cols]\n",
    "    x_test.loc[:, x_cols_to_scale] = x_scaler.transform(x_test.loc[:, x_cols_to_scale])\n",
    "    y_test = test_set.loc[:, y_cols].values\n",
    "\n",
    "    # Convert to PyTorch tensors\n",
    "    x_train_tensor = torch.tensor(x_train.values, dtype=torch.float32)\n",
    "    y_train_tensor = torch.tensor(y_train, dtype=torch.float32)\n",
    "    x_test_tensor = torch.tensor(x_test.values, dtype=torch.float32)\n",
    "\n",
    "    return x_train_tensor, y_train_tensor, x_test_tensor, y_test, y_scaler\n",
    "\n",
    "def learn(model, x_batch, y_batch, criterion, optimizer):\n",
    "    \"\"\"Computes the model estimate, compares it to the ground-truth \n",
    "    and performs the parameter update.\n",
    "    \"\"\"\n",
    "    optimizer.zero_grad() \n",
    "    output = model(x_batch)  # forward pass through the model\n",
    "    loss = criterion(output, y_batch)\n",
    "    loss.backward()  # backward step to compute gradients\n",
    "    optimizer.step()  # optimizer uses gradient for parameter update\n",
    "\n",
    "    return loss.item()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "input_dim = x_train.shape[1]\n",
    "criterion = nn.MSELoss()\n",
    "\n",
    "# Training parameters\n",
    "batch_size = 64\n",
    "epochs = 10\n",
    "\n",
    "y_cols = ['id_k1', 'iq_k1', 'i_norm_k1']\n",
    "x_cols = [c for c in ohe_df if c not in y_cols]\n",
    "x_cols_to_scale = [c for c in x_cols if not c.startswith(('n_k', 'n_1k'))]\n",
    "\n",
    "for i, (train_index, test_index) in enumerate(kf.split(ohe_df)):\n",
    "    # Initialize model and optimizer inside loop\n",
    "    model = FeedForwardNetwork(input_dim)\n",
    "    optimizer = optim.Adam(model.parameters(), lr=5e-4)\n",
    "\n",
    "    # Train-test split\n",
    "    train_set = ohe_df.loc[train_index, :]\n",
    "    test_set = ohe_df.loc[test_index, :]\n",
    "    \n",
    "    x_train_tensor, y_train_tensor, x_test_tensor, y_test, y_scaler = prepare_data(\n",
    "        train_set, test_set, y_cols, x_cols, x_cols_to_scale\n",
    "    )\n",
    "\n",
    "    # DataLoader\n",
    "    train_dataset = TensorDataset(x_train_tensor, y_train_tensor)\n",
    "    train_loader = DataLoader(train_dataset, batch_size=batch_size, shuffle=True)\n",
    "\n",
    "    # Training loop\n",
    "    model.train()\n",
    "    for epoch in range(epochs):\n",
    "        progress_bar = tqdm(train_loader, desc=f\"Epoch {epoch+1}/{epochs}\")\n",
    "        for x_batch, y_batch in progress_bar:\n",
    "            loss = learn(model, x_batch, y_batch, criterion, optimizer)\n",
    "            progress_bar.set_postfix(loss=loss)\n",
    "    \n",
    "    # Prediction\n",
    "    model.eval()\n",
    "    with torch.no_grad():\n",
    "        prediction = model(x_test_tensor)\n",
    "        prediction = y_scaler.inverse_transform(prediction.numpy())\n",
    "        print(f\"MSE: id_k1: {criterion(torch.tensor(prediction[:,0]), torch.tensor(y_test[:,0])):.2f} A²,\\\n",
    "              iq_k1: {criterion(torch.tensor(prediction[:,1]), torch.tensor(y_test[:,1])):.2f} A²,\\\n",
    "              i_norm_k1: {criterion(torch.tensor(prediction[:,2]), torch.tensor(y_test[:,2])):.2f} A²\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbgrader": {
     "grade": true,
     "grade_id": "cell-d5d9aba25ae5a099",
     "locked": false,
     "points": 0,
     "schema_version": 3,
     "solution": true,
     "task": false
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Residuals\n",
    "fig, axes = plt.subplots(nrows=1, ncols=3, figsize=(10, 3), sharey='row')\n",
    "error = y_test - prediction\n",
    "\n",
    "for i in range(len(y_cols)):\n",
    "    axes[i].scatter(y_test[:, i], error[:, i], s=.5, label=y_cols[i])\n",
    "    axes[i].grid(alpha=.5)\n",
    "    axes[i].axhline(ls='--')\n",
    "    if i == 0:\n",
    "        axes[i].set_ylabel('Error in A')\n",
    "    axes[i].set_xlabel(y_cols[i]+' in A')\n",
    "fig.tight_layout()"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Create Assignment",
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.19"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
